#############################################################################
# Summary table and figures
#
#############################################################################

rm(list=ls())
graphics.off()

library(data.table)
library(fixest)
library(xtable)
library(Hmisc)
library("RColorBrewer")
library(ggplot2)
library(ggpattern)

rm(list=ls())
save_res = "Results"
time0_fig = Sys.time()

#############################################################################
# 1. Bring data
#############################################################################

#######################
### FE
#######################

load("Data_new/FE_results.RData")

## Marginal effects per county-year
data_FE = copy(dt_res)
data_FE[,mg1:=NULL]
setnames(data_FE,c("mg2","mg3"),c("fe_hom","fe_het"))
summary(data_FE)
dim(data_FE)

rm(list=setdiff(ls(),list("dt_map","data_FE","time0_fig")))

#######################
### JK
#######################

load("Data_new/JK2_results.RData")
data_JK = copy(pdf_sim)
dim(data_JK)
summary(data_JK)

# Keep separately the bias-corrected mean and variance
jk_mean = J0$R["Mean","corr"]
jk_var = J0$R["Var","corr"]

rm(list=setdiff(ls(),list("dt_map","data_FE","data_JK","jk_mean","jk_var","time0_fig")))

#######################
### GFE
#######################

load("Data_new/GFE/GFE_results.RData")

# Keep relevant variables
names_now=c("fips","year","beta","group","mg")
dc_g2 = dc_g2[,.SD,.SDcols=c(names_now,"ddr","corn_area")]
dc_g3 = dc_g3[,.SD,.SDcols=names_now]
dc_g4 = dc_g4[,.SD,.SDcols=names_now]
dc_g5 = dc_g5[,.SD,.SDcols=names_now]
dc_g6 = dc_g6[,.SD,.SDcols=names_now]
dc_g7 = dc_g7[,.SD,.SDcols=names_now]
dc_g8 = dc_g8[,.SD,.SDcols=names_now]

# Put the results in one data
names_now=c("beta","group","mg")
setnames(dc_g2,names_now,paste0(names_now,"_g2"))
setnames(dc_g3,names_now,paste0(names_now,"_g3"))
setnames(dc_g4,names_now,paste0(names_now,"_g4"))
setnames(dc_g5,names_now,paste0(names_now,"_g5"))
setnames(dc_g6,names_now,paste0(names_now,"_g6"))
setnames(dc_g7,names_now,paste0(names_now,"_g7"))
setnames(dc_g8,names_now,paste0(names_now,"_g8"))
dc = merge(dc_g2,dc_g3,by=c("fips","year"),all=TRUE)
dc = merge(dc,dc_g4,by=c("fips","year"),all=TRUE)
dc = merge(dc,dc_g5,by=c("fips","year"),all=TRUE)
dc = merge(dc,dc_g6,by=c("fips","year"),all=TRUE)
dc = merge(dc,dc_g7,by=c("fips","year"),all=TRUE)
dc = merge(dc,dc_g8,by=c("fips","year"),all=TRUE)
dim(dc)
summary(dc)

# Data at county-year level of marginal effects
data_GFE =dc[,.SD,.SDcols=c("fips","year",paste0("mg_g",2:8))]
setkey(data_GFE,"fips","year")
data_GFE[,cor:=.I]
dim(data_GFE)
summary(data_GFE)

# Data at the county level of groups and coefficients
data_GFE_gr =dc[,.SD,.SDcols=c("fips","year",paste0("beta_g",2:8),paste0("group_g",2:8))]
data_GFE_gr =data_GFE_gr[,lapply(.SD,mean),by="fips"]
summary(data_GFE_gr)
dim(data_GFE_gr)

rm(list=setdiff(ls(),list("dt_map","data_FE","data_JK","data_GFE","data_GFE_gr","jk_mean","jk_var","time0_fig")))

#######################
### RE 
#######################

load("Data_new/RE_results.RData")

### Correlated RE
RE_1 = copy(res_w2$dt_mg_sim)
RE_1 = RE_1[,.(id,x,post_den)]
setnames(RE_1,c("x","post_den"),c("x_re_cor","den_mg_re_cor"))

### Distribution of posterior means
RE_2 = copy(res_w2$dt_ext)
setnames(RE_2,c("id","t"),c("fips","year"))
RE_2[,id:=.I+max(RE_1$id)]

### Put together 
data_RE = merge(RE_1,RE_2,by="id",all=TRUE)
dim(data_RE)
summary(data_RE)

res_cre = res_w2
rm(list=setdiff(ls(),list("dt_map","data_FE","data_JK","data_GFE","data_GFE_gr","data_RE","res_cre","jk_mean","jk_var","PM2","time0_fig")))

#######################
### FE per period
#######################

load("Data_new/FE_period_results.RData")

## Marginal effects per county-year
data_FE_pp = copy(dt_res)
data_FE_pp = data_FE_pp[,.SD,.SDcols=c("fips","period","year","mg_fep")]
summary(data_FE_pp)
dim(data_FE_pp)

rm(list=setdiff(ls(),list("dt_map","data_FE","data_JK","data_GFE","data_GFE_gr","data_RE","data_FE_pp","res_cre","jk_mean","jk_var","time0_fig")))

#######################
### GFE per period
#######################

load("Data_new/GFE_period/GFE_periods_results.RData")

# Keep relevant variables
names_now=c("fips","year","beta","group_new","mg")
dc_g2 = dc_G2[,.SD,.SDcols=c(names_now,"ddr","corn_area","period")]
dc_g3 = dc_G3[,.SD,.SDcols=names_now]
dc_g4 = dc_G4[,.SD,.SDcols=names_now]

# Put the results in one data
names_now=c("beta","group_new","mg")
names_new=c("beta","group","mg")
setnames(dc_g2,names_now,paste0(names_new,"_g2"))
setnames(dc_g3,names_now,paste0(names_new,"_g3"))
setnames(dc_g4,names_now,paste0(names_new,"_g4"))
dc = merge(dc_g2,dc_g3,by=c("fips","year"),all=TRUE)
dc = merge(dc,dc_g4,by=c("fips","year"),all=TRUE)
dim(dc)
summary(dc)

# Data at county-year level of marginal effects
data_GFE_pp =dc[,.SD,.SDcols=c("fips","year","period",paste0("mg_g",2:4))]
setkey(data_GFE_pp,"fips","year")
data_GFE_pp[,cor:=.I]
dim(data_GFE_pp)
summary(data_GFE_pp)

rm(list=setdiff(ls(),list("dt_map","data_FE","data_JK","data_GFE","data_GFE_gr","data_RE","data_FE_pp","data_GFE_pp","res_cre","jk_mean","jk_var","time0_fig")))

#######################
### Put together
#######################

data_all = merge(data_FE,data_GFE,by=c("fips","year"),all=TRUE)
data_RE[,cor:=.I+max(data_all$cor)]
data_all = merge(data_all,data_RE[,.(fips,year,w,mu1,x_re_cor,den_mg_re_cor,cor)],by="cor",all=TRUE)
summary(data_all[,fips.y-fips.x])
summary(data_all[,year.y-year.x])
setnames(data_all,"fips.x","fips")
setnames(data_all,"year.x","year")
data_all[,fips.y:=NULL]
data_all[,year.y:=NULL]

data_JK[,cor:=.I+max(data_all$cor)]
summary(data_JK)
data_all = merge(data_all,data_JK[,.(cor,cdf_mont)],by="cor",all=TRUE)
setnames(data_all,"cdf_mont","mg_jk")
dim(data_all)
summary(data_all)

data_all = merge(data_all,data_FE_pp,by=c("fips","year"),all=TRUE)

setnames(data_GFE_pp,paste0("mg_g",2:4),paste0("mg_pp_g",2:4))
data_GFE_pp[,period:=NULL]
data_GFE_pp[,cor:=NULL]
data_all = merge(data_all,data_GFE_pp,by=c("fips","year"),all=TRUE)

### Review
dim(data_all)
names(data_all)
summary(data_all[!is.na(fips)])
summary(data_all[!is.na(mu1)])
summary(data_all[!is.na(x_re_cor)])
summary(data_all[!is.na(mg_jk)])

# Bring map US
load("Data_ori/Data_us_map.RData")

source("Code/0. CODE Auxiliary.R")
rm(data_FE,data_JK,data_GFE,data_RE,data_FE_pp,data_GFE_pp)
time1_fig=Sys.time()
save.image("Data_new/Data_summary.RData")


#############################################################################
# 2. Table of marginal effects
#############################################################################

rm(list=ls())
load("Data_new/Data_summary.RData")
save_results="Results/"

# Some distributions are already simulated with weights (JK)
data_all[,w:=corn_area]
data_all[is.na(w),w:=1]
summary(data_all[w==1])
summary(data_all[w!=1])

########################################
# Main time invariant models
########################################

nn = c("fe_hom","fe_het","mg_jk","mg_g5")
ss=sum_table(data=data_all,name_vars=nn,weight_var="w")

## Add column with RE
ss = cbind(ss,re_unc=res_cre$table_mg[1:7,])
ss

## Save
l=dim(ss)[2]
dig=rep(3,l+1)
print(xtable(ss,type="latex",digits=dig,display=c("s",rep("f",l))),hline.after = NULL,
      file="Results/Summary_table_time_inv_nose.tex",include.rownames = TRUE,include.colnames=FALSE,
      sanitize.text.function = function(x){x},only.contents = TRUE)

########################################
# GFE several
########################################

nn = paste0("mg_g",2:8)
sum_table(data=data_all,name_vars=nn,save_name="Summary_table_GFE",weight_var = "w")

########################################
# GFE pp several
########################################

nn = paste0("mg_pp_g",2:4)
data_all[period==1,paste0(nn,"_p1"):=.SD,.SDcols=nn]
data_all[period==2,paste0(nn,"_p2"):=.SD,.SDcols=nn]
data_all[period==3,paste0(nn,"_p3"):=.SD,.SDcols=nn]
nn0 = paste0("_p",1:3)
nn1 = c(paste0("mg_pp_g2",nn0),paste0("mg_pp_g3",nn0),paste0("mg_pp_g4",nn0))
summary(data_all[,.SD,.SDcols=nn1])
sum_table(data=data_all,name_vars=nn1,save_name="Table_time_varying_GFE",weight_var = "w")


#############################################################################
# 3. Densities
#############################################################################

data_all[year>=1990,fe_sub:=fe_het]
data_all[year>=1990,fe_short:=mgs]

# Basic structure
g = ggplot(data=data_all) + theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.text = element_text(size=18),axis.title = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = 'right',legend.text = element_text(size=14)) +
  scale_color_discrete(name= "Groups") +
  xlab("Temperature effects") + ylab("Density") 
yl=c(0,0.16)
xl=c(-30,10)

########################################
# FE weighted and unweighted
########################################

g + geom_density(aes(x=fe_het,weight=corn_area),adjust=1.5,linetype="dashed") +
  geom_density(aes(x=fe_het),adjust=1.5,linetype="solid") +
  xlim(xl) + ylim(yl)
ggsave(paste0(save_results,"Density_fe_weights.pdf"),width=6.3,height=3.9,units="in",dpi=300)

########################################
# FE long vs short
########################################

# Density pf marginal effects with homogeneous and heterogeneous. Weighted
g + geom_density(aes(x=fe_het,weight=corn_area),adjust=1.5) + 
  geom_density(aes(x=fe_short,weight=corn_area),adjust=1.5,linetype="dashed") + 
  geom_vline(xintercept = mean(data_all$fe_hom,na.rm=TRUE),linetype="longdash") +
  xlim(c(-30,20)) + ylim(yl)
ggsave(paste0(save_results,"Density_fe_long_vs_short.pdf"),width=6.3,height=3.9,units="in",dpi=300)


########################################
# FE vs JK2
########################################

g + geom_density(aes(x=fe_het,weight=corn_area),adjust=1.5,linetype="dashed") +
  geom_density(aes(x=mg_jk),adjust=3) +
  xlim(xl) + ylim(yl)
ggsave(paste0(save_results,"Density_jk2.pdf"),width=6.3,height=3.9,units="in",dpi=300)

########################################
# RE 
########################################

g + geom_line(data=data_all,aes(x=x_re_cor,y=den_mg_re_cor)) + 
    geom_density(data=data_all,aes(x=fe_het,weight=corn_area),adjust=1.5,linetype="dashed") +
    xlim(xl) + ylim(yl)
ggsave(paste0(save_results,"Density_RE.pdf"),width=6.3,height=3.9,units="in",dpi=300)

########################################
# FE per period 
########################################

save_results = "Results/"

breaks=1:3
p_color = palette.colors(palette = "Tableau10")[1:7]
p_color = c(p_color[3],p_color[5],p_color[1])
p_color = setNames(p_color,breaks)
p_color

ggplot(data=data_all,aes(x=mg_fep,weight=corn_area,color=factor(period),linetype=factor(period),size=factor(period))) + theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.title = element_text(size=18),axis.text = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = c(0.2,0.70),legend.text = element_text(size=14),
        legend.key.size=unit(1.5,'line'),legend.key.width = unit(3, "line"), 
        legend.box.background = element_rect(colour = "black",size=1)) +
  xlab("Temperature effects") + ylab("Density") + 
  scale_color_manual(name = "Subperiod", breaks = c("1","2","3"),
                     labels = c("1950-1968","1969-1987","1988-2005"),
                     #values = c("1" = "#F8766D","2" = "#00BA38","3" = "#00A9FF")) +
                     values = p_color) +
  scale_linetype_manual(name = "Subperiod", breaks = c("1","2","3"),
                        labels = c("1950-1968","1969-1987","1988-2005"),
                        values = c("1" = "dashed","2" = "dotted", "3" = "dotdash")) +
  scale_size_manual(name = "Subperiod", breaks = c("1","2","3"),
                    labels = c("1950-1968","1969-1987","1988-2005"),
                    values = c("1" = 0.8,"2" = 0.8, "3" = 0.8)) +
  xlim(xl) + ylim(yl) +
  stat_density(geom="line",position="identity",adjust=2) +
  geom_density(aes(x=fe_het,color="",size="",linetype=""),adjust=1.5,linetype="solid",color="black",size=0.5) 
ggsave(paste0(save_results,"Density_FE_per_period.pdf"),width=6.3,height=3.9,units="in",dpi=300)

########################################
# GFE
########################################

dt_plot = data_all[is.na(corn_area)==FALSE,.SD,.SDcols=c("fips","year",paste0("mg_g",2:8),"corn_area")]
summary(dt_plot)
dim(dt_plot)
dt_plot = reshape(dt_plot,direction="long",idvar=c("fips","year"),timevar="G",varying=paste0("mg_g",2:8),sep="")
dt_plot = copy(dt_plot)
summary(dt_plot)
dim(dt_plot)
dim(dt_plot)[1]/7
dt_plot[,tot_area:=lapply(.SD,sum),by="G",.SDcols="corn_area"]
dt_plot[,w:=corn_area/tot_area]

p_breaks = c("2","3","4","5","6","7","8")
p_types = c("2"="stripe","3"="stripe","4"="none","5"="stripe","6"="stripe","7"="crosshatch","8"="circle")
p_angles = c("2"=0,"3"=45,"4"=135,"5"=0,"6"=45,"7"=135,"8"=0)
p_spacing = c("2"=0.01,"3"=0.01,"4"=0.01,"5"=0.02,"6"=0.02,"7"=0.02,"8"=0.01)
p_color = setNames(palette.colors(palette = "Tableau10")[1:7],p_breaks)

g = ggplot(data=dt_plot[G<=8],aes(x=mg_g,weight=w,color=factor(G),fill=factor(G),pattern=factor(G),pattern_spacing=factor(G),pattern_angle=factor(G))) + 
  geom_bar_pattern(alpha=0.4,position="dodge",pattern_density=0.1,width=0.5,linewidth=0.1)  + theme_bw() +
  scale_pattern_manual(name= "Groups",breaks=p_breaks,values=p_types) +
  scale_pattern_angle_manual(name= "Groups",breaks=p_breaks,values=p_angles) +
  scale_pattern_spacing_manual(name= "Groups",breaks=p_breaks,values=p_spacing) +
  scale_color_manual(name= "Groups",breaks=p_breaks,values = p_color) + 
  scale_fill_manual(name= "Groups",breaks=p_breaks,values = p_color) +
  xlab("Temperature effects") + ylab("") +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.text = element_text(size=18),axis.title = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = c(0.1,0.55),legend.text = element_text(size=14),
        legend.key.size=unit(1.5,'line'),legend.box.background = element_rect(colour = "black",size=1))
g + geom_density(data=data_all[is.na(fips)==FALSE],
                 aes(x=fe_het,weight=corn_area,color=NULL,fill=NULL,pattern=NULL,pattern_spacing=NULL,pattern_angle=NULL),alpha=0,color="black",adjust=1.5,linetype="solid") +
  xlim(xl) #+ ylim(yl)
ggsave(paste0(save_results,"Histogram_GFE.pdf"),width=6.3,height=3.9,units="in",dpi=300)

########################################
# GFE per period
########################################

### Histograms per period, several values of G (same sample as FE per period)
dt_plot = data_all[is.na(year)==FALSE,.SD,.SDcols=c("fips","year","period","corn_area",paste0("mg_pp_g",2:4))]
table(dt_plot[,.(year,period)])

nn = paste0("mg_pp_g",2:4)
dt_plot[period==1,paste0(nn,"_p1"):=.SD,.SDcols=nn]
dt_plot[period==2,paste0(nn,"_p2"):=.SD,.SDcols=nn]
dt_plot[period==3,paste0(nn,"_p3"):=.SD,.SDcols=nn]
summary(dt_plot)

dt_plot = reshape(dt_plot,direction="long",idvar=c("fips","year"),timevar="G",varying=paste0("mg_pp_g",2:4),sep="")
dt_plot = copy(dt_plot)
dim(dt_plot)
summary(dt_plot)
table(dt_plot$G)
dt_plot[,tot_area:=lapply(.SD,sum),.SD="corn_area",by=c("G","period")]
dt_plot[,w:=corn_area/tot_area]
yl=c(0,0.6)

#### Color blind color and patterns
p_breaks = c("2","3","4")
p_color = setNames(palette.colors(palette = "Tableau10")[c(6,3,1)],p_breaks)
p_types = c("2"="none","3"="circle","4"="stripe")
p_angles = c("2"=45,"3"=45,"4"=45)
p_spacing = c("2"=0.01,"3"=0.01,"4"=0.02)

# Period 1
ggplot(data=dt_plot[period==1],aes(x=mg_pp_g,weight=w,color=factor(G),fill=factor(G),pattern=factor(G),pattern_spacing=factor(G),pattern_angle=factor(G))) + 
  geom_bar_pattern(alpha=0.4,position="dodge",width = 0.5,linewidth=0.1)  + theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.text = element_text(size=18),axis.title = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = c(0.1,0.7),legend.text = element_text(size=14),
        legend.key.size=unit(1.5,'line'),legend.box.background = element_rect(colour = "black",size=1)) +
  scale_pattern_manual(name= "Groups",breaks=p_breaks,values=p_types) +
  scale_pattern_angle_manual(name= "Groups",breaks=p_breaks,values=p_angles) +
  scale_pattern_spacing_manual(name= "Groups",breaks=p_breaks,values=p_spacing) +
  scale_color_manual(name= "Groups",breaks=p_breaks,values = p_color) + 
  scale_fill_manual(name= "Groups",breaks=p_breaks,values = p_color) + 
  xlab("Temperature effects") + ylab("") +
  xlim(xl) + ylim(yl) 
ggsave(paste0(save_results,"Histogram_GFE_per_period_p1.pdf"),width=6.3,height=3.9,units="in",dpi=300)

# Period 2
ggplot(data=dt_plot[period==2],aes(x=mg_pp_g,weight=w,color=factor(G),fill=factor(G),pattern=factor(G),pattern_spacing=factor(G),pattern_angle=factor(G))) + 
  geom_bar_pattern(alpha=0.4,position="dodge",width = 0.5,linewidth=0.1)  + theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.text = element_text(size=18),axis.title = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = c(0.1,0.7),legend.text = element_text(size=14),
        legend.key.size=unit(1.5,'line'),legend.box.background = element_rect(colour = "black",size=1)) +
  scale_pattern_manual(name= "Groups",breaks=p_breaks,values=p_types) +
  scale_pattern_angle_manual(name= "Groups",breaks=p_breaks,values=p_angles) +
  scale_pattern_spacing_manual(name= "Groups",breaks=p_breaks,values=p_spacing) +
  scale_color_manual(name= "Groups",breaks=p_breaks,values = p_color) + 
  scale_fill_manual(name= "Groups",breaks=p_breaks,values = p_color) + 
  xlab("Temperature effects") + ylab("") +
  xlim(xl) + ylim(yl) 
ggsave(paste0(save_results,"Histogram_GFE_per_period_p2.pdf"),width=6.3,height=3.9,units="in",dpi=300)

# Period 3
ggplot(data=dt_plot[period==3],aes(x=mg_pp_g,weight=w,color=factor(G),fill=factor(G),pattern=factor(G),pattern_spacing=factor(G),pattern_angle=factor(G))) + 
  geom_bar_pattern(alpha=0.4,position="dodge",width = 0.5,linewidth=0.1)  + theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.text = element_text(size=18),axis.title = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = c(0.1,0.7),legend.text = element_text(size=14),
        legend.key.size=unit(1.5,'line'),legend.box.background = element_rect(colour = "black",size=1)) +
  scale_pattern_manual(name= "Groups",breaks=p_breaks,values=p_types) +
  scale_pattern_angle_manual(name= "Groups",breaks=p_breaks,values=p_angles) +
  scale_pattern_spacing_manual(name= "Groups",breaks=p_breaks,values=p_spacing) +
  scale_color_manual(name= "Groups",breaks=p_breaks,values = p_color) + 
  scale_fill_manual(name= "Groups",breaks=p_breaks,values = p_color) + 
  xlab("Temperature effects") + ylab("") +
  xlim(xl) + ylim(yl) 
ggsave(paste0(save_results,"Histogram_GFE_per_period_p3.pdf"),width=6.3,height=3.9,units="in",dpi=300)


#############################################################################
# 4. Maps
#############################################################################

rm(list=ls())
load("Data_new/Data_summary.RData")
source("Code/0. CODE Auxiliary.R")
save_res = "Results"

# Aggregate marginal effects at the county level (weighted) for those that are not weighted already 
dt_plot = data_all[is.na(fips)==FALSE,lapply(.SD,wtd.mean,weights=corn_area),by="fips",.SDcols=c("fe_hom","fe_het",paste0("mg_g",2:8),"mg_fep",paste0("mg_pp_g",2:4))]
dim(dt_plot)
summary(dt_plot)
qq= quantile(dt_plot$fe_het,probs=seq(0.05,0.95,0.05))
qq
bmin= qq[2]
bmax= qq[18]

# FE long panel
plot_map(data=dt_plot[,.(fips,fe_het)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_fe_het",palette="OrRd",rev=1,nc=9)

# GFE 
for (k in seq(2,8,1)) {
  plot_map(data=dt_plot[,.SD,.SDcols=c("fips",paste0("mg_g",k))],trunc=1,bmin=bmin,bmax=bmax,save_name=paste0("Map_mg_gfe_g",k),palette="OrRd",rev=1,nc=9)
}

# RE Posterior means 
plot_now = res_cre$dt[,.(id,mu1)]
setnames(plot_now,"id","fips")
plot_map(data=plot_now,trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_re_post_means",palette="OrRd",rev=1,nc=9)

# FE per period
dt_plotp = data_all[!is.na(year),lapply(.SD,wtd.mean,weights=corn_area),.SDcols=c("mg_fep",paste0("mg_pp_g",2:4)),by=c("fips","period")]
dim(dt_plotp)
summary(dt_plotp)
plot_map(data=dt_plotp[period==1,.(fips,mg_fep)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_period_p1",palette="OrRd",rev=1,nc=9)
plot_map(data=dt_plotp[period==2,.(fips,mg_fep)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_period_p2",palette="OrRd",rev=1,nc=9)
plot_map(data=dt_plotp[period==3,.(fips,mg_fep)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_period_p3",palette="OrRd",rev=1,nc=9)

# GFE per period G=4
plot_map(data=dt_plotp[period==1,.(fips,mg_pp_g4)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_G4_per1",palette="OrRd",rev=1,nc=9)
plot_map(data=dt_plotp[period==2,.(fips,mg_pp_g4)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_G4_per2",palette="OrRd",rev=1,nc=9)
plot_map(data=dt_plotp[period==3,.(fips,mg_pp_g4)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_G4_per3",palette="OrRd",rev=1,nc=9)

#############################################################################
# 5. Figures intro
#############################################################################


rm(list=setdiff(ls(),list("data_all","save_res","time0_fig")))
g= ggplot(data=data_all)  +  theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), axis.text = element_text(size=18),
        axis.title = element_text(size=18),
        legend.position = 'none',legend.text = element_text(size=14)) 

##################################
# Sparcity
##################################

### Histogram with a big mass point at zero and a few non-zeros ("lasso") 90%
set.seed(1)
dt_sim = data.table(id=seq(1,1000),v0=runif(n=1000)*3-2)
dt_sim[,v1:=v0]
dt_sim[id<=900,v1:=0]
summary(dt_sim)
dt_sim[abs(v1)<0.98,v1:=0]
dt_sim[v1>=0.98,v1:=v1/2]
w=6.3
h=3.9

g + geom_histogram(data=dt_sim,aes(x=v1,y=..density..)) +
  scale_x_continuous(breaks=0,lim=c(-2.2,1)) + xlab(expression(paste('Coefficient ', beta[i]))) +
  scale_y_continuous(breaks=0,lim=c(0,9)) +ylab("Histogram")
ggsave(paste0(save_res,"/Intro_1a.pdf"),width=w,height=h,units="in",dpi=300)

### A discrete histogram ("groups") 
vaux = sort(runif(n=10)*3-2)
dt_sim[,v2:=vaux[1]]
dt_sim[,dist:=abs(v0-v2)]
for (k in seq(2,10)) {
  dt_sim[,dist_new:=abs(v0-vaux[k])]
  dt_sim[dist_new<dist,v2:=vaux[k]]
  dt_sim[dist_new<dist,dist:=dist_new]
}
table(dt_sim$v2)
dt_sim[v2>0.6,v2:=1.2]

g + geom_histogram(data=dt_sim,aes(x=v2,y=..density..)) +  
  scale_x_continuous(breaks=0,lim=c(-2.2,1))  + xlab(expression(paste('Coefficient ', beta[i]))) +
  scale_y_continuous(breaks=0,lim=c(0,9))  +ylab("Histogram")
ggsave(paste0(save_res,"/Intro_1b.pdf"),width=w,height=h,units="in",dpi=300)

### A continuous normal density ("random effects") 
g + stat_function(data=dt_sim,aes(x=v0),fun=function(x) dnorm(x,mean=-0.9,sd=0.5)) +
  scale_x_continuous(breaks=0,lim=c(-2.2,1)) + xlab(expression(paste('Coefficient ', beta[i]))) +
  scale_y_continuous(breaks=0,lim=c(0,0.9))  +ylab("Density")
ggsave(paste0(save_res,"/Intro_1c.pdf"),width=w,height=h,units="in",dpi=300)

##################################
# Over time
##################################

Nsim = 30
TT=50
set.seed(1)
dt_sim=data.table(id=rep(seq(1,Nsim),TT),t=sort(rep(seq(0,TT-1),Nsim)))
dt_sim[t==1,beta:=runif(n=Nsim)]
dt_sim[,beta:=mean(beta,na.rm=TRUE)-1.05,by="id"]
dt_sim[id==1,delta:=-0.05*(t+1)]
dt_sim[,delta:=mean(delta,na.rm=TRUE),by="t"]
dt_sim[,v1:=beta+delta]

dt_sim[id==1,delta2:=(0.05*t+1)]
dt_sim[,delta2:=mean(delta2,na.rm=TRUE),by="t"]
dt_sim[,v2:=beta*delta2]

dt_sim[,group:=1*(id<=Nsim/3)+2*(id>Nsim/3 & id<=2*Nsim/3)+3*(id>2*Nsim/3)]
dt_sim[group==1, v3:= delta*0.2+min(dt_sim$beta)-2]
dt_sim[group==2, v3:= -1.2*delta*min(dt_sim$beta)-0.1]
dt_sim[group==3, v3:= delta*0.2+max(dt_sim$beta)]
summary(dt_sim)

ylim=c(-4,0)

# 2a: beta_i+delta_t 
g + geom_line(data=dt_sim,aes(x=t,y=v1,color=factor(id))) +
  xlab('Time t') + ylab(expression(paste('Coefficient ', beta[it]))) +
  scale_y_continuous(breaks=0,lim=ylim) + scale_x_continuous(breaks=NULL) +
  scale_color_manual(values=rep("black",length(unique(dt_sim$id))))
ggsave(paste0(save_res,"/Intro_2a.pdf"),width=w,height=h,units="in",dpi=300) 

# 2b: beta_i*delta_t 
g + geom_line(data=dt_sim,aes(x=t,y=v2,color=factor(id))) +
  xlab("Time t") + ylab(expression(paste('Coefficient ', beta[it]))) +
  scale_y_continuous(breaks=0,lim=ylim) + scale_x_continuous(breaks=NULL) +
  scale_color_manual(values=rep("black",length(unique(dt_sim$id))))
ggsave(paste0(save_res,"/Intro_2b.pdf"),width=w,height=h,units="in",dpi=300)

# 2c: beta_(k_i,t) 
g + geom_line(data=dt_sim,aes(x=t,y=v3,color=factor(id))) +
  xlab("Time t") + ylab(expression(paste('Coefficient ', beta[it]))) +
  scale_y_continuous(breaks=0,lim=ylim) + scale_x_continuous(breaks=NULL) +
  scale_color_manual(values=rep("black",length(unique(dt_sim$id))))
ggsave(paste0(save_res,"/Intro_2c.pdf"),width=w,height=h,units="in",dpi=300)

time1_fig = Sys.time()

